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ABSTRACT 

The dissolution process of star clusters is rather intricate for theory. We investigate it 
in the context of chaotic dynamics. We use the simple Plummer model for the gravi- 
tational field of a star cluster and treat the tidal field of the Galaxy within the tidal 
approximation. That is, a linear approximation of tidal forces from the Galaxy based 
on epicyclic theory in a rotating reference frame. The Poincare surfaces of section 
reveal the effect of a Coriolis asymmetry. The system is non-hyperbolic which has 
important consequences for the dynamics. We calculated the basins of escape with 
respect to the Lagrangian points Li and L2. The longest escape times have been mea- 
sured for initial conditions in the vicinity of the fractal basin boundaries. Furthermore, 
we computed the chaotic saddle for the system and its stable and unstable manifolds. 
The chaotic saddle is a fractal structure in phase space which has the form of a Cantor 
set and introduces chaos into the system. 
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1 INTRODUCTION 

The dissolution process of star clusters is an old problem 
in stellar dynamics. Once a star cluster has formed some- 
where in a galaxy, it tends to lose mass due to dynamical 
interactions until it has completely dissolved. It turned out 
that the physics behind the dissolution of star clusters is 
intricate and fascinating for theory. If a star cluster of fi- 
nite mass were isolated, in virial equilibrium (i.e. {v1) — 
12cri£,) and the velocity distribution given by a Maxwellian 
fM{X) = (4/y^) exp(-X2), where X = v/(V2aio) and 
v,Ve and am are the velocity, the escape velocity and the 
velocity dispersion, respectively, the fraction of stars which 
are faster than the rms escape speed were given by Xe = 
J^fM{X)dX = 2y67^exp(-6) + erfc^e ~ 0.00738316. 
This simple analytical result was published by Ambartsum- 
ian (1938) and two years later, independantly by Spitzer 
(1940) who named this effect "evaporation". The relevant 
process which brings stars above the escape speed and lets 
them evaporate, is two-body relaxation. The time scale of 
relaxation, which determines the rate of dynamical evolu- 
tion of a star cluster, yields thus an upper limit to the 
lifetime of any star cluster. However, since Xe is so small 
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(and relaxation time relatively long), the evaporation time 
is much longer than a Hubble time for typical globular clus- 
ters. Following a suggestion of Chandrasekhar (1942), King 
(1959) studied the effect of "potential escapers". These are 
stars which have been scattered above the escape energy 
but which have not yet left the cluster. These may be scat- 
tered back to negative energies within a crossing time and 
remain bound. Henon (1960, 1969) stressed the importance 
of few close encounters between stars for the rate of mass loss 
of star clusters. The Fokker-Planck approximation, which 
is widely used to study the dynamical evolution of star 
clusters, neglects strong encounters by construction. Nev- 
ertheless, close encounters could still be interpreted statis- 
tically as a certain discontinuous Markov process (Tschar- 
nuter 1971). However, direct A'^-body models seem to be ide- 
ally suited to study this phenomenon in more detail. Spitzer 
& Shapiro (1972) estimate that " close encounters may pro- 
duce effects perhaps as great as 10 percent of the "dominant" 
distant encounters" and ignore them. 

Nature provides an environment for star clusters in 
which the escape rate is typically strongly enhanced as com- 
pared with the slow evaporation rate of isolated star clusters: 
The tidal field of a galaxy induces saddle-like troughs in the 
walls of the potential well of a star cluster (cf. Figure 1). It 
therefore lowers the energy threshold in star clusters above 
which stars can escape from zero to a negative value (Wielen 
1972, 1974). Moreover, if we consider a star cluster in the 
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tidal field of a galaxy as a dynamical system, the tidal field 
can change the system's dynamics in a dramatical way as 
compared with an isolated system. 

In general, the escape process from star clusters in a 
tidal field proceeds in two stages: (1) Scattering of stars 
into the "escaping phase space" by two-body encounters 
and (2) leaJcage through openings in the equipotential su- 
faces around saddle points of the potential. The "escaping 
phase space" is defined as the subset of phase space, from 
which escape is possible. It is well understood that the time 
scale for a star to complete stage (1) scales with rela^xation 
time. On the other hand, the time scale for a star to com- 
plete stage (2) depends mainly on its energy (but also on its 
location in phase space as we will sec). When wc neglect the 
effect of two-body relaxation for the consideration of stage 
(2), the motion of a single star in the star cluster is deter- 
mined between times ti and t2 only by the smooth grav- 
itational potential in which the star moves. The potential 
itself is generated by the other stars in the star cluster dis- 
regarding their "grainyness" and by the superposed galactic 
gravitational field, which is due to the matter distribution 
of the galaxy. Within this framework, we will study stage 
(2) of the escape process in this paper. In this connection, 
the work of Fukushige & Heggie (2000) is of major interest. 
Their main result is an expression for the time scale of es- 
cape for a star in stage (2) which has just completed stage 
(1). The dependance of the escape process on two time scales 
which scale differently with the particle number A'' imposes a 
severe scaling problem for A''-body simulations. The scaling 
problem is of relevance since the it is on today's general- 
purpose hardware architectures not yet simply feasible to 
simulate the evolution of globular clusters with realistic par- 
ticle numbers of a few hundred thousands or even millions of 
stars by means of direct A''-body simulations. The result of 
Fukushige & Heggie has been applied in Baumgardt (2001) 
to solve the important scaling problem for the dissolution 
time of star clusters in the special case of circular cluster 
orbits. The obtained scaling law tdis oc ij!^*, where tdis and 
trh are the dissolution and half-mass relaocation times, re- 
spectively has been verified, e.g. in Spurzem et al. (2005). 

The problem of escape has also a long history in the 
context of the theories of dynamical systems and chaos. 
It is well-known for a long time, that certain Hamiltonian 
systems allow for escape of particles towards infinity. Such 
"open" Hamiltonian systems have been studied by Rod 
(1973), Churchill ot al. (1975), Contopoulos (1990), Con- 
topoulos & Kaufmann (1992), Siopis et al. (1997), Navarro 
& Henrard (2001) and Schneider, Tel & Neufeld (2002). The 
related chaotic scattering process, in which a particle ap- 
proaches a dynamical system from infinity, interacts with 
the system and leaves it, escaping to infinity, was investi- 
gated by many authors, as Eckhardt & Jung (1986), Jung 
(1987), Jung & Scholz (1987), Eckhardt (1987), Jung & Pott 
(1989), Bleher, Ott & Grebogi (1989) and Jung & Ziemniak 
(1992). Chaotic scattering in the restricted three-body prob- 
lem has been studied by Benet et al. (1997, 1999). Typically, 
the infinity acts as an attractor for an escaping particle, 
which may escape through different exits in the equipoten- 
tial surfaces. Thus it is possible to obtain basins of escape (or 
"exit" basins), similar to basins of attraction in dissipativc 
systems or the well-known Newton- Raphson fractals. Special 
types of basins of attraction (i.e. "riddled" or "intermingled" 



basins) have been explored by Ott et al. (1993) and Som- 
merer & Ott (1993, 1996). Basins of escape have been stud- 
ied by Bleher et al. (1988), and they are discussed in Con- 
topoulos (2002). Reasearch on escape from the paradigmatic 
Hcnon-Hciles system has been done by dc Moura & Letc- 
her (1999), Aguirre, Vallejo & Sanjuan (2001), Aguirre & 
Sanjuan (2003), Aguirre, Vallejo & Sanjuan (2003), Aguirre 
(2004) and Seoane Seoane Sepiilveda (2007). These papers 
served as the basis of this work. Relatively early, it was rec- 
ognized, that the key to the understanding of the the chaotic 
scattering process is a fractal structure in phase space which 
has the form of a Cantor set (Cantor 1884) and is called the 
chaotic saddle. Its skeleton consists of unstable periodic or- 
bits (of any period) which arc dense on the chaotic saddle 
(e.g. Lai 1997) and introduce chaos into the system (e.g. 
Contopoulos 2002). The properties of chaotic saddles have 
been investigated by different authors, as Hunt (1996), Lai 
et al. (1993), Lai (1997) or Mottcr & Lai (2001). Both hyper- 
bolic and non-hyperbolic chaotic saddles occur in dynamical 
systems. In the first case, there are no Kolmogorov-Aruold- 
Moser (KAM) tori, which means that all periodic orbits are 
unstable. In the second case, there are both KAM tori and 
chaotic sots in the phase space (J. C. Vallejo, priv. comm. 
and e.g. Lai et al. 1993). We note that all of the above ref- 
erences on the chaotic dynamics are exemplary rather than 
exhaustive since there exists a vast amount of literature on 
these topics. 

The aim of this paper is to allude to the importance of 
this last-mentioned branch of research for the field of stellar 
dynamics of star clusters. We will study the escape process 
from star clusters within the framework of chaotic dynam- 
ics. In section 2, we introduce the tidal approximation, i.e. 
approximate equations of motion for stellar orbits in a star 
cluster which is embedded in the tidal field of a galaxy. In 
section 3, we describe our (very simple) model of the gravi- 
tational potential. In section 4, wc discuss Poincarc surfaces 
of section, which show the effect of a Coriolis asymmetry. 
Furthermore, we discuss the basins of escape in Section 5 
and the chaotic saddle and its stable and unstable invariant 
manifolds in Section 6. Section 7 contains the discussion and 
conclusions. 



2 THE TIDAL APPROXIMATION 

The "tidal approximation" which is widely used in stellar 
and galactic dynamics for studies of stellar systems in a tidal 
field is nothing else than a simple approximation which, his- 
torically, has been applied already in the 19th century in 
"Hill's problem" (e.g. Stumpff 1965, Szebehely 1967, Siegel 
& Moser 1971) in the context of the (rather intricate) lu- 
nar theory. The difference to the tidal approximation lies 
merely in the form of the gravitational potentials which are 
used. The assumption that a star cluster moves around the 
Galactic centre on a circular orbit allows to use the epicyclic 
approximation to calculate steady linear tidal forces acting 
on the stars in the star cluster. As in the circular restricted 
three-body problem the appropriate coordinate system is a 
rotating reference frame in which both the star cluster cen- 
tre and the Galactic centre (i.e., the primaries) are at rest. 
Its origin is the star cluster centre, sitting in the minimum of 
the effective Galactic potential. The x-axis points away from 
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the Galactic centre; the y-axis points in the direction of the 
rotation of the star cluster around the Galactic centre; the 
z-axis lies perpendicular to the orbital plane and points to- 
wards the Galactic North pole. We define scaled corotating 
coordinates (x, y, z) as 



x = {R-Rg)/rt, y ^ Rg 



■jt) /rt, z = z jrt (1) 



where (i?, 0, z') are galactocentric cylindrical coordinates, 
Rg and are the radius and the frequency of the circular 
orbit, respectively, rt is a length scale (i.e. the tidal radius 
defined in Equation (111 below) and t is time (in the context 
of "Hill's problem" cf. Glaschke 2006). Since the coordinate 
system is rotating, centrifugal and Coriolis forces appear ac- 
cording to classical mechanics. In addition, tidal forces enter 
the equations of motion for stellar orbits near the origin of 
coordinates. To first order, we have in the rotating frame 



where 
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is the effective potential, which is plotted in Figure 1 for the 
2D case. Other isolating integrals are not given in the form 
of a simple analytical expression. However, some solutions 
of ([gJi - Q are subject to a third integral (Henon & Heiles 
1964), as has been demonstrated numerically by Fukushige 
& Heggie (2000), who calculated a Poincare surface of sec- 
tion. In principle, one may obtain power series expansions 
of such third integrals, see, e.g. the original works by Gus- 
tavson (1966) and Finkler, Jones & Sowell (1990) and the 
review in Moser (1968). Also, third integrals can be related 
to the existence of Killing tensor fields which are well-known 
in General Relativity (Clementi & Pettini 2002). At last, the 
tidal radius (King 1962) 
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where ^ci{x, y, z) and $g(-R, z') are the star cluster potential 
and the axisymmetric galactic potential, respectively. The 
second-last term on the right hand side in |2| is the cen- 
trifugal force and the last terms in (|2| and ^ are Coriolis 
forces. According to Binney & Tremaine (1987), the epicyclic 
frequency k and the vertical frequency v are given by 



' > (5) 



(flg,0) \ / (Hg,0) 

Thus the equations of motion can be written as 

X — fx — {n^ — 4uj'^)x + 2iuy (6) 
y = fy-2ux (7) 
z = fz- v'z, (8) 

where (/a;, Z^) = — V^ci is the (specific) force vector from 
the other cluster member stars which typically depends non- 
linearly on the coordinates. It is of interest for the following 
discussion that the equations of motion ([6]l - ([8| are invariant 
under time reversalQ Note that under a time reversal the 
frequencies also change their sign. Also, the equations of 
motion (|6| - ([8| admit an isolating integral of motion, the 
Jacobian 



C 



{i? 4- -I- i^) + $cff (a;, y, 



(9) 



^ The invariance under time reversal is related to the existence 
of a discrete group with only two elements, which acts on the 
space of solutions of the equations of motion l[6]l - (|8]l. In our 
case, the effective potential ( |10[ l and the Coriolis forces are time- 
symmetric, which implies the same symmetry of the equations of 
motion. 



rt 



(11) 



where Md is the star cluster mass, provides a fundamental 
length scale of the problem. It is the distance from the origin 
of coordinates to the Lagrangian points Li and L2 (which 
lie on the x axis, see Figure [T|. 



3 THE MODEL 

The characteristic frequencies u, k and v arise as proper- 
ties of the galactic gravitational potential. Throughout the 
paper, we use the values of the characteristic frequencies 
in the solar neighborhood. All of them can be expressed 
in terms of Oort's constants A and B (see, e.g., Binney 



Tremaine 1987): 



{A - BY 



= -AB{A - B), 



- 4lu^ = -4A{A - B), = A-KGpg + 2{A^ - B^). The 
vertical frequency v can be derived from the Poisson equa- 
tion for an axisymmetric system (see Oort 1965) and pg is 
the local Galactic density, which contributes to the dom- 
inant first term. We obtain both dimensionless parameters 
/uj^ ~ 1.8 and v'^ /uP ~ 7.6 using the values of Oort's con- 
stants given in Feast & Whitelock (1997) and the value for 
local Galactic density given in Holmberg & Flynn (2000). It 
is then convenient to choose the following system of units: 



G = 1, 



Mrf = 2.2 



(12) 



The resulting length unit is the tidal radius rt and the for- 
mulation of the dynamical problem with its equations of 
motion is completely dimensionless. For the star cluster, 
we use a Plummer model, the most simple analytic model 
for a star cluster (see appendix A). The density profiles of 
King models which have a cutoff radius, where the density 
drops to zero, fit the measured density profiles of globu- 
lar clusters better than Plummer models (King 1966). Since 
they are tidally limited by construction, they are at first 
glance ideally suited for our purpose. However, the grav- 
itational force field can only be tabulated from a numer- 
ical integration of a non-linear differential equation. We 
have made a compromise which is not relevant for the in- 
teresting physics: We choose the Plummer radius in such 
a way, that the Plummer model (see Appendix A) is the 
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Figure 1. Effective potential in the tidal approximation (2 = plane). The Lagrangian points at Li = (—1,0) and L2 = (1,0) can be 
seen. The escapers pass these saddle points while they leak out. The equipotential lines connecting them mark the tidal boundary of the 
star cluster. The details of the model are given in Section 3. 



best fit to a King model with Wo = 4 (i.e. with concen- 
tration c = logjo(rt/rK) — 0.840), which completely fills 
the Roche lobe in the tidal field, i.e. the density of the King 
model approaches zero at the tidal radius ( |11[ ). The fit of the 
density profiles is quite acceptable for density contrasts of 
logj^Q(pc/p('')) ^ 3, where pc and p{r) are the central density 
and the density as a function of radius, respectively. For the 
deviation between the King density profile and the Plummer 
fit we obtain {ppi{r) — pK{r))/pc < 1.2%. The ratio of the 
Plummer radius to the King radius and the "concentration" 
of the Plummer model (which can only be defined because 
of the existence of a tidal radius) are then rpi/vK — 1.257 
and cpi = logj^g{rt/rpi) ~ 0.741, respectively. In our units, 
the Plummer radius is therefore rpi ~ 0.182. 

As a more technical remark, we note that we used an 
8th-order Runge-Kutta scheme for the orbit integrations. 
The relative error in the Jacobian C was always limited to 
AC/C < 10"^^ for all orbit integrations. 



POINCARE SURFACES OF SECTION AND 
THE CORIOLIS ASYMMETRY 



A critical Jacobi constant 



Cl = 'l>cff(rt,0,0) = 'l'eff(-rt,0,0) 



(13) 
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Figure 3. The two main types of orbits at C = 0. Left: Regular 
retrograde orbit. Right: Chaotic prograde orbit. The variable C 
is defined in Equation l |14| . 



is given by the value of the effective potential ( 10 1 at the 
Lagrangian points Li and L2. For our model we have Cl ~ 
—3.264444506. For a Jacobian C > Cl the equipotential 
surfaces are open and particles can escape. Furthermore, we 
define the dimensionless deviation from Cl by 



C = (Cl - C)/Cl 



(14) 



where C is some other value of the Jacobian. The dimen- 
sionless deviation C is positive for C > Cl if C and Cl 
are both negative, which is always the case in this paper. 
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Figure 2. Poincare surfaces of section. Top left: At C = for orbits crossing y = witii y > 0, Top righit: Same as top left, but at 
C = 0.1, Bottom left: At C = fo r or bits crossing x = with y > 0, Bottom right: At C = for orbits crossing y = with i > 0. The 
variable C is defined in Equation (|l4|. 



A first insight can be gained by calculating Poincare sur- 
faces of section which are shown in Figure [2] for two dif- 
ferent Jacobi constants: The upper left surface of section 
is at the critical Jacobian Cl at which all orbits still re- 
main within a bounded area in phase space and we have 
no escapers. We can see that this is a system with divided 
phase space, i.e. we have both chaotic and regular orbits. It 
is striking that the left half of the surface of section is al- 
most completely occupied by regular, quasiperiodic orbits. 
These orbits are retrograde with respect to the orbit of the 
star cluster around the galactic centre (Fukushige & Heggie 
2000), as can be seen by looking at the sketch in Figure [4j 
and they are subject to a third integral of motion. On the 
other hand, most of the prograde orbits are chaotic apart 
from a few smaller regular islands. Since the effective poten- 



tial is mirror-symmetrical with respect to both the x- and 
y-axes (see Figure [T]), the asymmetry seen in the surfaces 
of section must be due to the Coriolis forces. Thus such a 
behaviour might be termed a "Coriolis asymmetry" (cf. In- 
nanen 1980). The Coriolis forces are special in the sense that 
their direction is not perpendicular to the tangent plane to 
the equipotential surfaces but to the velocity of a particle. 
The upper right surface of section is at a higher Jacobi con- 
stant. The particles can leak out through the openings in the 
equipotential surfaces and escape towards infinity (positive 
a;-direction) or the galactic centre (negative x-direction) . It 
is remarkable, that only the chaotic orbits escape, while the 
regular, quasiperiodic orbits remain within the tidal bound- 
ary of the star cluster, since the third integral restricts their 
accessible phase space and hinders their escape. In star clus- 
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Figure 4. Sketch of the coordinate system. The escapers leak out 
through the openings in the equipotential surfaces passing either 
L\ or L2. Only schematically, two orbits are shown which cross 
the X axis with y > 0. 



ters, two-body relaxation may scatter stellar orbits beyond 
the critical Jacobi constant. However, if the orbits are retro- 
grade, the stars will remain bound to the star cluster with 
high probability until two-body relaxation further scatters 
them into the escaping phase space. The two lower surfaces 
of section in Figure [2] indicate the orbital structure in po- 
sition space at C = Cl, where all orbits are restricted to 
the region within the almond-shaped tidal boundary of the 
star cluster (cf. Figure [T]). One notes that certain parts of 
the chaotic regions of the lower surfaces of section are less 
densely filled by stellar orbits which is a common property 
of dynamical systems. 

Figure [3] shows typical examples of the two main types 
of orbits. The regular orbit resembles a rosette orbit in an 
axisymmetric potential or a loop orbit which suggests that 
the third integral is some sort of generalization of angular 
momentum (Binney & Tremaine 1987). We found indeed 
that the angular momentum is approximately conserved for 
the left orbit of Figure [3] while it is not at all conserved for 
the right orbit. Other types of orbits can be found which 
are associated with the smaller regular islands in the sur- 
faces of section and which may look rather interesting. At 
last, we remark that the difference between retrograde and 
prograde orbits appears prominently in A'^-body models of 
dissolving rotating star clusters within the framework of the 
tidal approximation: Where either the regular or the chaotic 
domains in phase space are more strongly occupied by stel- 
lar orbits due to the existence of a net angular momentum 
of the cluster in one direction (see Ernst et al. 2007). 



5 THE BASINS OF ESCAPE 

Figure [5] shows the basins of escape for a tidally limited star 
cluster within the framework of the tidal approximation. 
For the plots on the left-hand side, ~ 3 x 10'' orbits have 
been integrated; for the plots on the right-hand side, it were 
~ 6 — 8 X 10^ orbits, depending on the area of the surface of 
section containing initial conditions. The phase space is di- 
vided into the escaping phase space (red and yellow regions) 
and the non-escaping phase space (black regions). The red 
regions denote initial conditions, where the escaping stars 
pass L\ while the yellow regions represent initial conditions, 
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Table 1. Fraction of orbits (in percent) belonging to the inter- 
section of the basins of escape with Poincare surfaces of section 
which are shown in Figure |5] 

where the escaping stars pass 1/2- The black regions show 
those initial conditions, where stars do not escape. These 
are first and foremost the regular regions where a third inte- 
gral is present. Note that the stable manifold of the chaotic 
saddle is not marked in black, since it is of Lebesgue measure 
zero (cf. Section 6). Since there exist regular regions (with 
KAM tori), the system is non-hyperbolic, i.e. there exist sta- 
ble periodic orbits with corresponding elliptic points in the 
surfaces of section (cf. Section 4 and Figure [2|. We can see 
that in the escaping phase space there exist regions, where 
we have a very sensitive dependance of the escape process 
on the initial conditions, i.e. a slight change of the initial 
conditions makes the star escape through the opposite La- 
grangian point. This is the classical indication of chaos. It 
is interesting to note that these regions arise from imme- 
diate vicinity of the black regions where orbits are regular. 
In these domains of phase space the red and yellow regions 
are completely intertwined with respect to each other: The 
boundary between these regions is fractal. The volume in 
phase space occupied by these regions (with sensitive depen- 
dance on the initial conditions) increases as the Jacobian C 
approaches the critical Jacobi constant Cl (i.e. in the limit 
C ^ 0) and the exits become smaller. At C ^ there is 
a maximal "fractalization" of phase space; when the sys- 
tem approaches the limit of small exits the basins become 
uncertain (Aguirre & Sanjuan 2003). The term "uncertain" 
means that we become unable to follow the real trajectory 
of a particle by means of numerical integrationj^ Moreover, 
the following theorem can be formulated: "_For all points P 
in the escaping phase space of an open Hamiltonian system, 
and for all S > (precision of the experiment), there exists 
a critical size of the exits Wc > such that for all w ^ Wc 
we can find a point P' in a hall centreed m P and radius S 
that belongs to a different basin than P" (Aguirre & Sanjuan 
2003). 

Table[l]shows the fraction of orbits (in percent) belong- 
ing to the intersection of the basins of escape with Poincare 
surfaces of section which are shown in Figure [5] It can be 
seen that in the limit (7^0 the fractions of particles pass- 
ing Li and L2 tend to be equal while this must not be the 
case if there are large areas without sensitive dependance of 
the escape process on the initial conditions. 

Figure [6] shows how the escape times are distributed 
on surfaces of section. The longest escape times correspond 
to initial conditions near the boundaries between the basins 

^ In other words, the computer fails here to be a Laplacian demon 
(Laplace 1814). 
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Figure 5. The basins of escape. Top left: At C = 0.1 for orbits crossing y = with y > 0, Top right: At C = 0.1 for orbits crossing 
y = with X > 0, Middle row: As the upper row, but at C = 0.01, Botton row: As the middle row, but at C = 0.001. The red regions 
denote initial conditions, where the escaping stars pass Li while the yellow regions represent initial conditions, where the e scap ing stars 
pass L2. The black regions show those initial conditions, where stars do not escape. The variable C is defined in Equation jl4|. 
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Figure 6. Distribution of escape times te on surfaces of section 
for C = 0.01. Top: For tiie x — Vx surface of section of Figurejs] 
Bottom: For tiie x — y surface of section of Figure |5] The darker 
tlie color, the longer the escape time. 
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Figure 7. Histogram of the fraction of remaining (non-escaped) 
orbits Ne{te > t)/A^e,o after time t. Top: For the x — surfaces 
of section of Figure [6] Bottom: For the x — y surfaces of section of 
Figure [e] Solid: 6 = 0.1, Dashed: C = 0.01, Dotted: C = 0.001. 
The inlays with two logarithmic axes show the early phase for the 
solid line (i.e. for C = 0.1). The variable C is defined in Equation 



long escape times which correspond to the regions with sen- 
sitive dependance on the initial conditions). The slopes of 
the exponentials (i.e. the decay constants) depend on the 
value of C but are identical for both surfaces of section. The 
exponential decay law indicates that the underlying process 
is of a statistical nature similar to the radioactive decay of 
unstable nuclides or that of bubbles in beer foam. 



of escape of Figure [5] The shortest escape times have been 
measured for the ordered regions without sensitive depen- 
dance on the initial conditions, i.e. those far away from the 
fractal basin boundaries. 

Figure [t] shows the fraction of remaining (non-escaped) 
orbits Ne{te. > t)/ Ne.fi after time t corresponding to the 
basins of escape shown in Figure [5] Only the escaping orbits 
have been used for the statistics. The orbits corresponding to 
the regions without sensitive dependance on the initial con- 
ditions shown in Figures [5] have short escape times as can be 
seen in Figure [6] For these orbits, the decay law is a power 
law as shown in the inlays for the solid curve (C = 0.1). On 
the other hand, the decay law is exponential for the chaotic 
orbits near the fractal basin boundaries (i.e. the orbits with 



6 THE CHAOTIC SADDLE AND ITS 
INVARIANT MANIFOLDS 

The stable manifold of the chaotic saddle is shown in the 
top row of Figure [8] for two Poincare surfaces of section. 
The stable manifold coincides with the fractal basin bound- 
aries of Figure [5] and therefore acts as a separatrix between 
the exit basins.. With data points of finite size, the top row 
shows orbits which do not escape for time i — > oo, although, 
strictly speaking, their Lebesgue measure is zero. The un- 
stable manifold of the chaotic saddle is shown in the middle 
row of Figure m These are orbits which do not escape for 
time t ~oo. Note that the stable and unstable manifolds 
are symmetric with respect to each other, since the equa- 
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tions of motion ([6| - ([8| are time-symmetric. For the plots 
in the middle row of Figure [S] the sign of the time step in 
the Runge-Kutta integrator has been reversed, The bottom 
row of Figure [8] shows the intersection of the chaotic sad- 
dle (i.e. a non-hyperbolic chaotic invariant set) with the two 
Poincare surfaces of section. It is the invariant set of non- 
escaping orbits for time t ~* oo and t — oo. The chaotic 
saddle has the form of a Cantor set (Cantor 1884) which is 
formed by the intersection of its stable and unstable man- 
ifolds. The fact that the system is non-hyperbolic implies 
that there are tangencies between the stable and unstable 
manifolds, i.e. that their angle is not always bounded away 
from zero (Lai et al. 1993). The unstable (hyperbolic) points 
in the intersection of the Poincare surface of section with the 
chaotic saddle correspond to unstable periodic orbits. As is 
well-known (e.g. Contopoulos 2002), these introduce chaos 
into the system since they repel the orbits in their neigh- 
bourhood in the direction of their unstable eigenvectors. 
The remarkable similarity of our system with the Henon- 
Heiles system (see Aguirre, Vallejo & Sanjuan 2001) is that 
the fractal dimension of the chaotic saddle tends to three 
(i.e. the black areas in Figure [s] grow until we have a maxi- 
mal fractalization of the phase space, cf. Aguirre & Sanjuan 
2003) , which is the dimension of the hypersurface of phase 
space with constant Jacobian. At (7 = there is a sudden 
transition where the non-hyperbolic invariant set abruptly 
fills the whole non-regular subset of phase space within the 
last closed equipotential surface and no escape is possible 
any more. This situation can be seen in the Poincare sur- 
faces of section in Figure [2] in Section 4 which shows the 
chaotic domains of phase space as a dotted area. 



7 DISCUSSION AND CONCLUSIONS 

We studied the chaotic dynamics within a star cluster which 
is embedded in the tidal field of a galaxy. We calculated 
within the framework of the tidal approximation Poincare 
surfaces of section, the basins of escape, the chaotic saddle 
and its stable and unstable manifolds as well. The system 
is non-hyperbolic which has important consequences for the 
dynamics, i.e. there are orbits which do not escape if relax- 
ation is neglected. These are mainly the retrograde orbits 
as has been shown earlier by Fukushige & Heggie (2000) 
and, more recently, in the iV-body study by Ernst et al. 
(2007). The escape times are longest for initial conditions 
near the fractal basin boundaries. The decay law is a power 
law for those stars which escape from the regions without 
sensitive dependance on the initial conditions in Figure 5] 
(i.e. with short escape times, as can be seen in Figure |6 1. 
On the other hand, the decay law is exponential for orbits 
which escape from the regions with sensitive dependance on 
the initial conditions (i.e. with long escape times). The ef- 
fect of relaxation (i.e. a diffusion in the Jacobian and the 
third integral among different stellar orbits) on the chaotic 
dynamics which we investigated in this work may be a very 
interesting topic for future research. 



8 ACKNOWLEDGEMENTS 

We thank Dr. Juan Carlos Vallejo and Prof. Burkhard Fuchs 
for pointing out some references to us, which were necessary 
for this work. Also, we thank Prof. Rudolf Dvorak for his 
comments on the manuscript. AE would like to thank his 
father for a discussion and Dr. Patrick Glaschke for many 
inspiring discussions. Furthermore, he gratefully acknowl- 
edges support by the International Max Planck Research 
School (IMPRS) for Astronomy and Cosmic Physics at the 
University of Heidelberg. 



REFERENCES 

Ambartsumian V. A., 1938, Ann. Leningrad State Univ., 
22, 19 

Aguirre J., Vallejo J. C, Sanjuan M. A. F., 2001, Phys. 
Rev. E, 64, 066208 

Aguirre J., Sanjuan M. A. F., 2003, Phys. Rev. E, 67, 
056201 

Aguirre J., Vallejo J. C, Sanjuan M. A. F., 2003, Int. J. 
Mod. Phys. B, 17, 4171 

Aguirre J., 2004, PhD thesis, Universidad Rey Juan Carlos, 
Spain 

Baumgardt H., 2001, MNRAS, 325, 1323 

Benet L., Trautmann D., Sehgman T. H., 1997, Cel. Mech. 

Dyn. Astron., 66, 203 

Benet L., Sehgman T. H., Trautmann D., 1999, Cel. Mech. 
Dyn. Astron., 71, 167 

Binney J., Tremaine S., 1987, Galactic Dynamics, Prince- 
ton Univ. Press 

Bleher S., Grebogi C, Ott, E., Brown R., 1988, Phys. Rev. 
A, 38, 930 

Bleher S., Ott E., Grebogi C, 1989, Phys. Rev. Lett., 63, 
919 

Cantor G, 1884, Acta Mathematica, 4, 381. English transl. 
in ed. Edgar G. A., 1993, Classics on Fractals, Addison- 
Wesley 

Chandrasekhar S., 1942, Principles of Stellar Dynamics, 
Univ. Chicago Press, Chicago 

ChurchiU R. C, PeceUi G., Rod D. L., 1975, J. Diff. Eq., 
17, 329 

Clementi C, Pettini M., 2002, Cel. Mech. Dyn. Astron., 
84, 263 

Contopoulos G, 1990, Astron. Astrophys., 231, 41 
Contopoulos G, Kaufmann D., 1992, Astron. Astrophys., 
253, 379 

Contopoulos G, 2002, Order and Chaos m Dynamical As- 
tronomy, Berlin (Springer) 

De Moura A. P. S., Letelier P. S., 1999, Phys. Lett. A, 256, 
362 

Eckhardt B., Jung C, 1986, J. Phys. A: Math. Gen., 19, 
L829 

Eckhardt B., 1987, J. Phys. A: Math. Gen., 20, 5971 
Ernst A., Glaschke P., Fiestas J., Just A., Spurzem R., 
2007, MNRAS, 377, 465 

Feast M., Whitelock P., 1997, MNRAS, 291, 683 

Finkler P., Jones C. E., Sowell G. A., 1990, Phys. Rev. A, 

42, 1931 

Fukushige T., Heggie D. C, 2000, MNRAS, 318, 753 
Glaschke P., 2006, PhD thesis, University of Heidelberg, 
http://www.ub.uni-heidelberg.de/archiv/6553 



Escape from the vicinity of fractal basin boundaries of a star cluster 11 



Gustavson F., 1966, AJ, 71, 670 

Henon M., 1960, Ann. Astrophys., 23, 668 

Henon M, Heiles C, 1964, AJ, 69, 73 

Henon M., 1969, A&A, 2, 151 

Holmberg J., Flynn C, 2000, MNRAS, 313,209 

Hunt B. R., Ott E., Yorke J. A., 1996, Phys. Rev. E, 54, 

4819 

Innanen K. A., 1980, Ap. J., 85, 81 

Jung. C, 1987, J. Phys. A: Math. Gen., 20, 1719 

Jung C, Scholz H. J., 1987, J. Phys. A: Math. Gen., 20, 

3607 

Jung C, Pott S., 1989, J. Phys. A: Math. Gen., 22, 2925 
Jung C, Ziemniak E., 1992, J. Phys. A: Math. Gen., 25, 
3929 

King I. R., 1959, AJ, 64, 351 
King I. R., 1962, AJ, 67, 471 
King I. R., 1966, AJ, 71, 64 

Lai Y.-C, Grebogi C, Yorke J. A., Kan I., 1993, Nonhn- 
earity, 6, 779 

Lai Y.-C, 1997, Phys. Rev. E, 56, 6531 
Laplace, 1814, Theorie analytique des probabihtes. In: 
Oeuvres Completes de Laplace, Volume VH., Paris, 1820 
( Gauthier- Villars) 

Moser J. K., 1968, Lectures on Hamiltonian Systems, Mem. 
AMS, 81, 1. Reprinted in eds. MacKay R. S. and Meiss J. 
D., 1987, Hamiltonian Dynamical Systems, Bristol (Adam 
Hilger) 

Motter A. E. Lai Y., 2001, Phys. Rev. E, 65, 015205 
Navarro J. F., Henrard J., 2001, A&A, 369, 1112 
Oort J. H., 1965, Stellar Dynamics, in: Galactic Structure, 
eds. A. Blaauw & M. Schmidt, Univ. Chicago Press, p. 455 
Ott E., Sommerer J. C., Alexander J. C. Kan I., Yorke J. 
A., 1993, Phys. Rev. Lett., 71, 4134 
Plummer H. C., 1911, MNRAS, 71, 460 
Rod D. L., 1973, J. DifT. Eq., 14, 129 

Schneider J., Tel T., Neufeld Z., 2002, Phys. Rev. E, 66, 
066218 

Seoane Sepulveda J. M., 2007, PhD thesis, Universidad Rey 
Juan Carlos, Spain 

Siegel C. L., Moser J. K., 1971, Lectures on Celestial Me- 
chanics, Springer Edition, Berlin, Heidelberg, New York 
Siopis C, Kandrup H. E., Contopoulos G., 1997, Cel. Mech. 
Dyn. Astron., 65, 57 

Sommerer J. C, Ott E., 1993, Nature, 365, 138 
Sommerer J. C, Ott E., 1996, Phys. Lett. A, 214, 243 
Spitzer L. Jr., Shapiro S. L., 1972, Ap. J., 173, 529 
Spurzem R., Giersz M., Takahashi K., Ernst A., 2005, MN- 
RAS, 364, 948 

Stumpff K., 1965, Himmelsmechanik, Vol. 2, VEB 
Deutscher Verlag der Wissenschaften, Berlin 
Szebehely V., 1967, Theory of orbits. The restricted prob- 
lem of three bodies. New York (Academic Press) 
Tscharnuter W. M., 1971, Ap&SS, 14, 11 
Wielen R., 1972, in Lecar M., ed., Cravitational N-body 
problem. Proceedings of lAU Colloquium 10, Reidel Pub- 
lishing Company, Dordrecht 

Wielen R., 1974, in Mavridis L. N., ed.. Stars and the 
Milky Way System. Proceedings of the First European As- 
tronomocal Meeting, Athens, Sept. 4-9, 1972., Vol. 2, p. 
326 



APPENDIX A: THE ISOLATED PLUMMER 
MODEL 

The most well-known self-consistent model of a spherically 
symmetric star cluster is the Plummer model (Plummer 
1911). It is given by 



$(j/) = -$o 



= — $0 COS I 



(1-1-1/2) 

P'^y) = PO , „.5/2 = PO COS^ 

with y — r/rpi ^ 0, (f) — arctan(j/) and the Plummer radius 



(Al) 
(A2) 
(A3) 



rpi = 



GM ( 3A4" 



$0 V^Trpo 



1/3 



(A4) 



where M is the total mass, — "l>o is the central potential and 
po the central density, all of them being finite. From | |A2[ ) 
we find the half mass radius 

yn = tan [arcsin (s"^''^)] ^ 1.30476603 (A5) 

The velocity dispersion can be obtained by integrating the 
Jeans equation of hydrostatic equilibrium. It is given by 



2. X 1 *o 

(T (j/) — , — cos ( 



6 yrTF 6 



Therefore the relation 



vl(y) = I2a\y), 



(A6) 



(A7) 



where = ^2|$(y)| is the escape speed from an isolated 
Plummer model, strictly holds at any radius (cf. Section 1 
of this paper) . One can show that the Plummer model is the 
only model which has this property. 

This paper has been typeset from a TjnjX/ M^rjX file prepared 
by the author. 



